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Abstract 

The diversity of fungi along environmental gradients has been little explored in contrast to plants and animals. 
Consequently, environmental factors influencing the composition of fungal assemblages are poorly understood. The aim of 
this study was to determine whether the diversity and composition of leaf and root-associated fungal assemblages vary 
with elevation and to investigate potential explanatory variables. High-throughput sequencing of the Internal Transcribed 
Spacer 1 region was used to explore fungal assemblages along three elevation gradients, located in French mountainous 
regions. Beech forest was selected as a study system to minimise the host effect. The variation in species richness and 
specific composition was investigated for ascomycetes and basidiomycetes assemblages with a particular focus on root- 
associated ectomycorrhizal fungi. The richness of fungal communities associated with leaves or roots did not significantly 
relate to any of the tested environmental drivers, i.e. elevation, mean temperature, precipitation or edaphic variables such as 
soil pH or the ratio carbon:nitrogen. Nevertheless, the ascomycete species richness peaked at mid-temperature, illustrating a 
mid-domain effect model. We found that leaf and root-associated fungal assemblages did not follow similar patterns of 
composition with elevation. While the composition of the leaf-associated fungal assemblage correlated primarily with the 
mean annual temperature, the composition of root-associated fungal assemblage was explained equally by soil pH and by 
temperature. The ectomycorrhizal composition was also related to these variables. Our results therefore suggest that above 
and below-ground fungal assemblages are not controlled by the same main environmental variables. This may be due to 
the larger amplitude of climatic variables in the tree foliage compared to the soil environment. 
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Introduction 

Forest microbial assemblages hold major roles in ecosystem 
functioning. However, the distribution patterns of fungal assem- 
blages are poorly understood because few studies have been 
performed at large geographical scales [1-2]. The climatic factors 
influencing the microbial richness and composition are equally still 
poorly understood compared to macroorganisms. The diversity of 
macroorganisms decreases with increased latitude [3] and, 
depending on the group, a hump-shaped distribution or a 
decrease in species richness with elevation is observed for plants, 
vertebrates and invertebrates [4-5]. These large-scale distribution 
patterns are of special interest in the context of climate change. 
Indeed, there is an increasing amount of evidence showing 
changes in plant communities subject to global warming, in 
particular along elevation gradients with a shift in the distribution 
of plants species notably more pronounced at higher elevations [6- 
7]- 



While it appears that there is a considerable fungal diversity, 
communities of fungi have been less studied compared with those 
of macroorganisms. This is due to the difficulty in describing 
microbial communities adequately. However, advances in molec- 
ular techniques, such as the recent high-throughput sequence- 
based technologies now allow far easier characterisation of fungal 
assemblages and improved estimation of fungal species richness. 
Today, sequence-based identification is recognized as a powerful 
method that has significantly improved our perception of fungi in 
a variety of environmental conditions and habitats [8-10], 
particularly along environmental gradients [11-13]. 

Microorganisms might not follow the elevational diversity 
patterns generally observed for macroorganisms. Recent findings 
suggest that bacteria diversity may not decrease with elevation 
[14-15], and may be higher at mid-elevation (850 m) [18]. Recent 
studies have shown no change of fungal richness to be associated 
with elevation [17-18]. However, a majority of studies have 
conclusively noted a decrease in ectomycorrhizal (EcM) diversity 
associated with elevation [12], [19-21]. In addition to the 
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diversity, the composition of fungal assemblages varies with 
elevation, as reported for fungal assemblages of the beech 
phyllosphere [11], mycorrhizal fungi [12], [19], [22] and fungal 
wood decomposers [17]. Climatic variables may explain part of 
these variations in diversity or composition of the fungal 
assemblages. Indeed, the mean annual temperature and precip- 
itation explained the observed patterns of EcM fungal richness and 
assemblage structure along elevation and latitudinal gradients 
[10], [12], [16]. However, other factors may drive these patterns, 
particularly in the case of EcM fungi [10], [16], [23]. Noteworthy, 
the host plant and the soil pH could be the major drivers of below- 
ground fungal community diversity and composition [12], [16], 
[19], [24], [25]. 

In this study, we aimed to determine whether the richness of 
root-associated and leaf fungal assemblage matches the elevation 
diversity gradient observed for the majority of organisms i.e. shows 
a decrease with elevation. We hypothesized that i) fungal species 
richness decreases with higher elevation and lower temperature 
and ii) leaf-associated fungal assemblages are more constrained by 
climatic variables than root-associated assemblages because the 
below-ground assemblage does not experience environmental 
changes as large as the above-ground assemblage does. 

For this purpose, richness and composition of the fungal 
assemblages were analyzed along elevation gradients in three 
different mountainous regions of France. We focused on beech 
dominated forests to alleviate a potential host effect. The leaf "or 
phyllosphere" fungi we studied were all species inhabiting both the 
surface and the interior of the leaves [9] , the root-associated fungi 
were those associated with fine roots. 

Materials and Methods 

Ethics Statement 

Field sampling was permitted by the Office National des Forets 
and the private owner of the St Nicolas forest site. No protected 
species was sampled. 

Site description and sampling design 

Three gradients were located in the Vosges, the Alps and the 
Pyrenees mountains in France. Forest stands were chosen at three 
different elevations in the Vosges and the Alps and at five 
elevations in the Pyrenees. We selected stands with mature and 
dominant (>50% of the basal area) beech trees (Fagus sybatica L.). 

For each elevation site, three plots of five contiguous beech trees 
were selected. The trees were approximately 5 meters apart within 
a plot with distances between each plot of approximately 
25 meters. Three leaves per tree, from different branches, were 
sampled between the 29 th June and the 9 th July 2010 as described 
in [11]. At the same time, three soil cores were sampled per tree 
(approximately 50 cm 3 with 10 cm depth) at 1 m from the trunk 
in the north, south-west and south-east directions. 

The leaves were placed in individual plastic bags containing 
10 ml of silica gel each (Sigma- Aldrich, St. Louis, MO, USA) to 
ensure that they completely dried within a few hours. The plastic 
bags were stored at 16°C before DNA extraction. Four discs per 
leaf were cut and were ground into a homogenous powder. For a 
given tree, the three soil cores were pooled and stored at 4°C 
before processing (maximum 20 days, without visible mould 
development). The entire root system was sorted from the soil with 
a 5-mm-mesh sieve and frozen until processing. The root systems 
were then gently washed in tap water to remove soil particles along 
with broken root tips that had been washed out. To select for EcM 
fungi, the thicker roots (over 3 mm diameter) were discarded and 
the fine-roots (including residual rhizospheric soil) were pooled, 



lyophilized and ground 10 sec in a ball mill. Soil analyses were 
carried out by the Laboratoire d' Analyse des Sols d'Arras (http:// 
www5.lille.inra.fr/las) for pH (water method), total C (NF ISO 
10694), total N (NF ISO 13878) and P contents [26]. The mean 
annual temperature and the annual precipitation at each elevation 
site were both obtained from the AURELHY method (Analysis 
Using the Relief for Hydrometeorology) for the 1971-2001 period 
[27] (Table 1). 

DNA isolation and 454-pyrosequencing 

The leaf samples were processed as previously described [11]. 
The total DNA was isolated with a CTAB phenol/chloroform/ 
isoamyl alcohol protocol. Three leaf DNA isolations from the 
three different leaves were performed per tree. For homogenized 
root samples, the total genomic DNA was isolated from 50 mg of 
fine root powder per sample using the DNeasy Plant Mini Kit 
(Qiagen, Courtaboeuf, France). For each root DNA sample, 
approximately 15 mg of polyvinylpolypyrrolidone (PVPP) was 
added to neutralize putative PCR-inhibitors. Manufacturers' 
instructions were then followed. 

The Internal Transcribed Spacer 1 (ITS1) region was targeted 
as a recognized DNA barcode for fungal identification [28]. The 
fungal ITS1 region was amplified from the metagenomic DNA 
with the forward ITS IF (CTTGGTCATTTAGAGGAAGTAA) 
and reverse ITS2 (GCTGCGTTCTTCATCGATGC) primers. 
Before amplification, the primers were adjoined by a self- 
generated multiplex identifier (MID) and the adaptor A 
(C GT ATC GC CTCCCTCGCGC CATC AG) or by only the 
adaptor B (CTATGCGCCTTGCCAGCCCGCTCAG) for 454 
Titanium pyrosequencing of multiplexed amplicons. The adaptor 
A and the MID were fused with ITS2 for leaf samples and with 
ITS IF for root samples. Sequencing was unidirectional and started 
with the A adaptor. 

The amplification of the rDNA ITS1 for the trees from the 
same plot was done with the same tag-encoded primer set. PCR 
mix preparation was done in 20 \xL. For leaf samples, the protocol 
is described in [1 1] and is similar to the one followed by root 
samples. The final concentrations were: 0.4 x enzyme buffer, 
1.5 mM MgC12, 1.12 mg/mL bovine serum albumin, 0.2 mil 
dNTP each, 0.4 mM each tagged primer (ITS IF and ITS2), 
0.05 U/uL Taq polymerase (Sigma-Aldrich, St Louis, MO, USA). 
DNA extracts from root samples were diluted 100 times. The 
GeneAmp system 9700 was used with 50% ramp to decrease the 
rate of temperature change during the cycles. The program was as 
follow: 5 min at 95°C, 35 cycles of 30 sec at 95°C, 40 sec at 53°C 
and 45 sec at 72°C, and a final step 7 min at 72°C. The PCR 
reactions were done for each sample separately ( 1 1 sites x 3 plots 
x5 trees) and amplicons from the same plot were then pooled 
[29]. Each PCR product was purified with the AMPure XP 
purification kit (AgenCourt Bioscience, Fullerton, CA, USA) for 
leaf samples or on a MultiScreen HTS PCR Plate MSNU030 
(Millipore Corporation, Billerica, MA, USA) for root samples. The 
latter procedure consisted of two washings with 100 uL H 2 0 and 
final resuspension in 40 U.L H 2 0. The concentration of each 
sample was measured with the NanoDrop ND-1000 spectropho- 
tometer and the samples were pooled in equal proportion. In total, 
we obtained independent libraries corresponding to 66 leaf and 
root tagged samples, sequenced with three independent runs. The 
454 pyrosequencing was done at the Genoscope (Evry, France). 
The raw data were deposited on the Sequence Read Archive 
website (http://www.ncbi.nlm.nih.gov/sra) under the study ac- 
cession ERP003510. The three libraries are under the following 
experiment accessions: ERX280868 (alias TCA.AQN_GOTS_ 
G0YLLAW02) for the root samples, ERX280870 (alias 
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TCA.AQN_FOTS_GYUGVSB04) for the leaf samples from the 
Pyrenees, and ERX28087 1 (alias TCA.AQN_HOTS_G2Q8WDE03) 
for the leaf samples from the Alps and the Vosges. 

Bioinformatics 

The sequences were sorted into different files according to their 
multiplex identifier (MID) using sfffile software, with no mis- 
matches allowed for the MIDs, which were then removed. The 
sequences were trimmed with trim.flows in Mothur version 1.26.0 
with the default parameters [30]. 

The ITS1 was extracted using Fungal ITS extractor version 2 
[31]. After extraction of the ITS1 region, the sequences were 
filtered based on a minimal length of 100 bp. Out of a total of 
370 449 sequences, 228 104 (61.6%), were retained for down- 
stream analyses. The Molecular Operational Taxonomic Units 
(MOTUs) were generated by clustering the sequences with Uclust 
version 3.0 (Usearch version 6.0.152) [32]. The similarity 
threshold was set at 97% as a suggested consensus in the literature 
[33]. The taxonomic assignation of a consensus sequence 
generated by Uclust was done against a curated GenBank 
database (July 2012) selecting for ribosomal sequences of 
eukaryotes and against unidentified sequences in www.ncbi.nlm. 
nih.gov using the Basic Local Alignment Search Tool (BLAST) 
algorithm version 2.2.23 [34], with the filter turned off (i.e. 
segments of the query sequence that have low compositional 
complexity were not masked off). This setting decreased the 
number of sequences with no BLAST hits. Ten BLAST hits were 
considered when available. The phylum name was extracted from 
the NCBI taxonomy corresponding to the Genlnfo identifier (GI 
number). An additional BLASTN was carried out online for 
MOTUs with no hit identified. This allowed us to determine the 
phylum for few MOTUs. For the assignation at the genus level, 
only BLAST results for which identity >90% and e-value <e 
were considered. The percentage of identity, given in the 
BLASTN result, was the ratio of the length of identical bases 
over the length of the alignment. To increase the robustness of the 
assignation, the MOTUs were assigned to the genus if at least 
eight blast hits were congruent with each other for the genus name 
or if all would be congruent when fewer hits were available. This 
threshold was chosen arbitrarily but allowed us to avoid selecting 
the best blast hit which may have been misleading. The species 
name was determined only for MOTUs with a genus assignation 
and for which the percentage identity of a BLAST hit was equal or 
greater than the ITS1 variability among phylum. The ITS1 
variability among phylum was calculated from the Table S3 
provided by Nilsson and collaborators (2008) [35]. The ITS1 
variability for each phylum was calculated as the mean ITS 1 intra- 
specific variability within each phylum weighted by the number of 
sequences (29 502 sequences in all in [35]) and was as follows: 
Ascomycota, 2.88%; Basidiomycota, 4.98%; Chytridiomycota, 
7.81%; Glomeromycota, 9.48%. The mean ITS1 intra-specific 
variability within the Zygomycota group was 4.13%. Some 
synonyms were corrected by referencing the website Species 
Fungorum (http://www.speciesfungorum.org). The MOTUs were 
assigned to the species level if at least six BLAST hits were 
congruent with each other for the species name or all should be 
congruent when fewer hits were available. The threshold was 
again chosen arbitrarily. The MOTUs taxonomically assigned to a 
g E same species or matching to the same GI number at >90% 
identity and for an e-value <e~ 45 , were subsequently grouped 
| £ together reducing further the total number of MOTUs. Many 
MOTUs were then discarded according to the following criteria: 
1) non-fungal MOTUs; 2) MOTUs supported by only one 
sequence (singletons) or 3) MOTUs with no hit found in the 
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additional BLASTN and supported by less than 10 sequences. The 
two last criteria may be indicative of possible sequencing artefacts, 
and this is why these MOTUs were not used in the analysis. 

Statistical analyses 

The leaf and the root-associated assemblages were analysed 
separately. An assemblage was defined as the addition of all fungal 
molecular operational taxonomic units (MOTUs, a reasonable 
proxy for species) living in a particular niche (roots or leaves). 
Ascomycete and basidiomycete assemblages of both roots and 
leaves were also analysed separately, as well as a functional group, 
the root-associated EcM fungi. Whether a particular species was a 
root-associated EcM fungi was based on current literature [36- 
38]. The list of the MOTUs is given in Table SI. The table 
indicates the phylum and if the MOTUs is an EcM fungi. 

Richness was used as diversity measure and was defined by the 
estimated number of MOTUs per sample. It was estimated, by 
rarefaction, at the lowest sequencing depths available, which were 
1400 and 500 sequences per sample for the phyllosphere and the 
root-associated assemblages respectively. The rarefaction curves 
were computed for each sample (i.e. each plot and each habitat) 
with the rarefy function of the vegan R package [39]. Before 
rarefaction, one sequence per MOTU was randomly discarded to 
exclude the singletons from the dataset [42] .To estimate the 
richness of a particular group i, such as the ascomycetes, the 
basidiomycetes, or the EcM fungi, we made the hypothesis that the 
ratio of MOTUs belonging to the group i within the unknown 
MOTUs (i.e. those for which the phylum was unknown, 10%, 
Table 2) was the same as within the assigned MOTUs. The 
estimated richness for the group i in sample j (R 5 ) was thus 
calculated as follows: 

1 v assigned 

Ny, the number of MOTUs assigned to group i in sample j 

•^assigned! me number of MOTUs assigned at least to the phylum 
level in sample j 

ERj, the estimated fungal richness of the sample j, estimated at 
the sequencing depth mentioned above. 

The correspondence analysis (CoA) was used to discriminate the 
samples based on their species composition. The abundance 
matrix was rarefied at sequencing depth mentioned above with the 
rrarejy function of the vegan R package and then transformed into 
a presence/absence matrix. From this presence /absence matrix, 
the MOTUs detected in less than three plots were discarded to 
prevent distortion of the CoA by the rare MOTUs. The CoA was 
done with the dudi.coa function of the ade4 R package with default 
parameters [40]. 

The relationship between environmental variables was studied 
with a principal component analysis (PCA) and a Pearson's 
correlation table was used to select the variables that were not 
correlated for further analysis. The PCA was performed with the 
PCA function of the FactoMineR library. The Pearson's correla- 
tion coefficient and the p-value of the tests were obtained with the 
cor.test function of the stats library. To test for significant effects of 
selected environmental variables (mean annual temperature and 
soil pH) on fungal richnesses, mixed linear models were used with 
the Ime function available from the nlme R package [41]. The 
interactions were not included in the models because they were 
not significant. To test for significant effects of the selected 
environmental variables on specific compositions, permutational 
analyses of variance (PERMANOVA) were performed with the 



adonis function available from the vegan R package. For each 
model, the elevation site was nested within the region factor (Alps, 
Vosges or Pyrenees). The permutations were restricted within the 
region using the strata option of the adonis function. The 
significance level was set at 0.05 for each test. 

Results 

Description of the dataset 

We obtained 228 104 sequences that passed the quality filters. 
The mean number of sequences per sample following the cleaning 
steps was 4 544 and 2 368 for the leaf and root samples, 
respectively. 

There were 6 475 MOTUs generated by the clustering step 
based on the sequence alignment. The second clustering based on 
the taxonomic assignation led to 4 855 MOTUs. The MOTUs 
supported by only one sequence within the entire dataset 
accounted for a high percentage (37.8%) and will be hereafter 
referred to as singletons. The number of singletons per sample was 
significantly correlated to the number of MOTUs per sample 
(F,, 64 = 157.50, p-value<0.001, r 2 = 0.707). The singletons were 
discarded, leaving 3 018 non-singleton MOTUs (hereafter MO- 
TUs) for the analysis (Table 2). 

The MOTUs were found in one (49.3%), two (28.8%) or three 
(2 1 .9%) studied regions. Information on how many MOTUs were 
shared between and within regions is given in Table S2 and Figure 
SI). 29.0% were found in one out of the 11 sites and 1.5% were 
found in all the 1 1 elevation sites. A significant positive correlation 
was noted between the logarithm of the number of sequences per 
MOTUs and the frequency i.e. the number of plots in which the 
MOTUs were present (F 13016 = 5721.0, p-value<0.001, 
r 2 = 0.655). 

The MOTUs were unassigned (10.1%), assigned at the phylum 
(57.2%), the genus (15.3%) or the species (17.4%) level only. 
Altogether, Dikarya represented 86.4% of the MOTUs and 96.1% 
of the phylum-assigned MOTUs (Table 2). 

Taxonomic diversity of leaf and root-associated fungal 
assemblages 

We found a total of 1 457 nonsingletons MOTUs in the leaf 
samples and 1 701 MOTUs in the root samples (Table 2). These 
MOTUs represented a high diversity of fungal genera as 165 and 
2 1 2 different fungal genera were detected in the phyllosphere and 
root-associated assemblages, respectively. The most diverse genera 
were Cryptococcus (25 MOTUs) and Taphrina (22 MOTUs) in the 
phyllosphere, and Russula (55 MOTUs) and Cortinarius (50 
MOTUs) in the root samples. The phyllosphere harboured a 
higher proportion of ascomycetes compared to the root system 
(t 64 = 18.96, p-value<0.001). There were 140 MOTUs shared by 
the leaf and root habitats but they were not equivalendy 
represented in both habitats. In total, 306 MOTUs were classified 
as EcM fungi (Table 2) and belonged to 39 different genera. 
Among the EcM MOTUs, 289 (94.4%) were basidiomycetes. 

Environmental variations along the elevation gradients 

According to the correlation table (Table S3), the environmental 
variables (mean annual temperature, annual rainfall, soil carbon 
content, soil nitrogen content, carbon:nitrogen ratio, and soil 
phosphorus content) were all found to significantly correlate to 
elevation with the exception of soil pH. Moreover, with the 
exception of soil pH and C:N ratio, all environmental variables 
were significantly related to each other. The soil pH was found to 
be only weakly related to the C:N ratio. The C:N ratio was also 
significandy related to the soil phosphorus content. The PCA 
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Table 2. Leaf and root samples description. 







Combined dataset 


Leaf samples 


Root samples 




(66 samples) 


(33 samples*) 


(33 samples*) 


Total number of 454 reads 


370 449 


248 299 


122 150 


Number of quality sequences 


228 104 


149 946 


78 158 


Range of sequences per sample 




1 506-15 830 


552-3 622 


Range of sequences per elevation site 




5 691-30 380 


4 716-10 040 


Number of MOTUs 1 


4 855 


2 224 


2 771 


Number of non singletons 2 


3 018 


1 457 


1 701 


Range of non singletons per sample 




1 34-487 


1 34-342 


Range of non singletons per site 




255-634 


375-568 


Number of nonsingletons assigned at the following taxonomic level 


Phylum 


1 725 


967 


813 


Genus 


525 


143 


415 


Species 


462 


168 


343 


Unknown 


306 


179 


130 


Number of nonsingletons assigned to the following phyla 


Ascomycota 


1 672 (55.4%) 


958 (65.8%) 


799 (47%) 


Basidiomycota 


935 (31.0%) 


320 (22.0%) 


667 (39.2%) 


Zygomycota 


97 (3.2%) 


0 


97 (5.7%) 


Other or incertain fungal lineages 


8 (0.3%) 


0 


8 (0.5%) 


Unknown 


306 (10.1%) 


179 (12.3%) 


1 30 (7.6%) 


Ecological mode 


EcM 3 


306 


0 


306 


Other or unknown 


2 712 


1 457 


1 395 



* Three plots xll elevation sites. 

Molecular Operational Taxonomic Unit (97% similarity); 

2 MOTUs supported by two sequences at least within the combined dataset; 

including putative EcM fungi. 

doi:1 0.1 371 /journal.pone.01 00668.t002 



confirmed these results. The two first axes of the PCA on 
environmental variables explained 49.3% and 20.2% of the 
variability, respectively (Figure S2). Nitrogen, phosphorus, tem- 
perature and rainfall were the main contributors to the first axis 
with approximately 20% each, whereas C:N ratio and pH 
contributed to the second axis for 55% and 34% respectively. 
Therefore, only two variables, the mean annual temperature and 
the soil pH, were tested for explaining fungal assemblage diversity 
and composition (Table 3 and Table 4). 

Variations in fungal richness along the elevation 
gradients 

For leaf samples, the estimated richness ranged from 103 to 199 
MOTUs per plot (81 to 154 ascomycetes and 12 to 55 
basidiomycetes). For root samples, the estimated richness ranged 
from 88 to 157 MOTUs per plot (36 to 80 ascomycetes, 31 to 66 
basidiomycetes, and 13 to 35 EcM MOTUs) (Table S4). 

The richness was not significandy related to elevation (— 0.78< 
t<2.08; 0.076<p-value<0.638; results not shown). The only 
exception was for the root ascomycete richness, which was 
significandy correlated to elevation (t=— 3.47, p-value = 0.010) 
(Figure S3). However, there was a mid-elevational peak in species 
richness for root-associated fungi in the Pyrenees and in the 
Vosges (Figure S3). 

No relationship was found between richness and mean annual 
temperature except for the root ascomycete richness (t = 3.28, p- 



value = 0.014) (Figure 1; Table 3). Additionally, regardless of the 
subgroup considered (ascomycetes, basidiomycetes, EcM fungi), 
no significant relationship was found linking richness and other 
environmental variables (soil pH, C, N, C:N, P, rainfall; data 
shown only for pH; Table 3). Similar conclusions were obtained 
when the diversity of fungal community was measured using the 
Shannon diversity index (data not shown). 

The mountainous region (introduced as a random factor in the 
statistical analyses) was a notably different source of variability 
depending on the fungal assemblage. While leaf-associated fungal 
richness varied significantly depending on the region (45.5-68.7% 
of the total variance, Table 3, likelihood ratio of 6.6-7.8 with 
0.005<p-value<0.060), root assemblage richness did not (likeli- 
hood ratio of about 10 8 , p-value over 0.999, Table 3). For the 
global root-associated fungal community, most of the variation in 
richness was induced by the random site factor (likelihood ratio of 
6.9—11.4 with p-value < 0.0 10). The within-site variation between 
the 3 plots of 5 trees represented between 28 and 44% of the total 
variation depending on the assemblage considered (Table 3). 

Variations in fungal composition along the elevation 
gradients 

The CoA showed that relatively littie variation in the species 
composition existed between the 3 plots within a site (Figures 2A 
and 2B). Still, the within-site variability can be measured by the 
residuals of the permutational multivariate analysis of variance 
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and represented about 28-30% for leaf-associated communities 
and slightly more, 37-39% for root-associated communities 
(Table 4). 

The phyllosphere fungal assemblage of the Pyrenees was very 
different from those of the Alps and the Vosges, which were quite 
similar to each other (Figure 2A). The region effect represented 
about 25% of the variability in the PERMANOVA (Table 4). The 
sites were distributed on the first and second axes according to the 
elevation, in particular for the Pyrenees, with low elevation sites 
having high values on the first axis and low values on the second 
axis (Figure 2 A). This distribution indicates that the assemblage in 
the high elevation sites may be different from the low elevation 
sites. The ten leaf MOTUs that contributed the most on the first 
axis of the CoA contributed for over 0.86% each. The 
permutational multivariate analysis of variance confirmed that 
temperature was the main environmental factor affecting the 
composition of leaf-associated fungal communities (Table 4, see 
also [11]). Mean site temperature represented about 15% of the 
variation of the leaf-associated communities compared to about 
5% for soil pH. 

Concerning the root assemblage, the regional effect accounted 
for 15% of the variability (Table 4). On Figure 2B which relates to 
root-associated fungi, the three plots from the two sites in the 
Pyrenees which had strong positive values on the second axis were 
the two lowest elevation sites from the Pyrenees region (Lourdes 
and Laveyron, Table 1). The permutational multivariate analysis 
of variance showed that soil pH and site mean temperature 
explained similar levels of the root-associated communities 
variability, whether the entire community was considered or only 
part of it (ascomycetes, basidiomycetes or EcM species). Site 
temperature explained about 8% of the variability while soil pH 
explained about 9% (Table 4). Some EcM species were only 
detected in the warmer sites (&12.0°C) e.g. Russula Solaris and 
Russula ingwa. On the contrary, Melanogaster sp2 and Clavulina 
amethystina were only detected in the colder sites (S7.5°C). Twenty- 
three Russula spp, 16 Cortinarius spp, 14 Hypochrea spp, 10 
Cryptococcus spp., eigth Lactarius spp., five Boletus spp., and other 
species from other genera were not found in sites with a soil pH>5. 

The EcM fungi belonged to 39 different genera. From 18 to 28 
were detected in a single site. Twenty-three of the observed EcM 
genera were detected in the three regions, 1 0 were in 2 regions and 
6 only in one region. Nine genera were detected in all of the 1 1 
elevation sites. 



Discussion 

Fungal richness does not decline along elevation 
gradients 

Our results show that no clear relationship exists between fungal 
richness and elevation across a broad range (from 131 to 1533 m) 
(Figure S3). This finding adds to studies that have reported no 
patterns or inconclusive patterns of fungal richness along elevation 
gradients [17-18], [22] but differs to previous studies which found 
lower richness levels in high-elevation sites compared to low or 
mid-elevation sites [12], [19]. A possible explanation for these 
differences is the relative homogeneity of host factor. For instance, 
the diversity of EcM hosts is known to actively influence the 
diversity of EcM assemblage [25], [43]. This factor was controlled 
along the three elevation gradients in our study (dominated beech 
forest), whereas it was not controlled in several previous studies 
[12], [19]. A poor estimate of the diversity of microbial community 
could also impede the detection of a clear pattern of richness along 
the gradients. However, similar conclusions are obtained using 
diversity index such as the Shannon index suggested in the 
literature [44]. 

In our study, the fungal richness does not correlate with the 
climatic and soil variables that covaried with elevation (mean 
annual temperature, annual rainfall, soil carbon content, soil 
nitrogen content, carbomnitrogen ratio, and soil phosphorus 
content) (Table 3). Our results are not in agreement with previous 
studies reporting that the mean annual temperature and precip- 
itation drive EcM species richness along elevation gradients [12], 
[16]. Interestingly, only ascomycete species richness peaked at 
mid-temperature (Figure ID). This richness distribution is 
significandy validated by a quadratic curve illustrating a mid- 
domain effect model (data not shown). Previously, different 
analyses of elevational diversity revealed patterns of mid-eleva- 
tional peaks in mammal species richness [45-46], illustrating this 
mid-domain effect [47]. Because this trend is not observed for 
basidiomycetes, it is possible that the richness of ascomycetes and 
the richness of basidiomycetes are not related to the same 
environmental variables. This finding suggests that considering the 
whole fungal assemblage might blur the relationship of sub- 
assemblage richness with environmental gradients. A confirmation 
of this finding is needed through additional studies. That is why we 
recommend studying the ascomycetes and the basidiomycetes 
independentiy. 



Table 3. Relationships between potential factors affecting leaf and root-associated fungal richness. 







Fixed variables 




Random factors* 








Mean annual temperature 


Soil pH 


Region 


Site 


Residual 


Leaf MOTUs 


-1.56 (0.162) 


-1.17 (0.257) 


25.43 (68.7) 


5.08 (2.7) 


16.39 (28.5) 


Leaf Ascomycetes 


-0.44 (0.674) 


-0.61 (0.647) 


1 2.40 (45.5) 


6.09 (11.0) 


12.12 (43.5) 


Leaf Basidiomycetes 


-1.88 (0.089) 


-0.92 (0.366) 


1 2.45 (66.7) 


3.22 (4.5) 


8.19 (28.9) 


Root MOTUs 


1.10 (0.309) 


-0.64 (0.528) 


0.003 (0) 


14.78 (68.9) 


9.94 (31.1) 


Root Ascomycetes 


3.28 (0.014) 


0.61 (0.547) 


2.11 (6.6) 


5.92 (51.7) 


5.32 (41.7) 


Root Basidiomycetes 


-0.32 (0.757) 


0.01 (0.991) 


0.001 (0) 


7.78 (55.3) 


7.00 (44.7) 


Root EcM 


0.20 (0.845) 


-0.38 (0.706) 


0.0006 (0) 


5.16 (59.0) 


4.30 (41.0) 



The T-values (p-values) of the linear mixed model are presented for the variables tested and the standard deviations (% of the total variance explained [random-h 
residual]) are presented for the random factors. Significant p-values in bold at 5%. The richness was estimated at 1 400 and 500 sequences per leaf and root samples 
respectively. 

* Standard deviations associated with the region and site nested within region as random factors (% of the total random plus residual variance explained). 
doi:1 0.1 371 /journal.pone.01 00668.t003 
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Table 4. Relationships between potential factors affecting leaf and root-associated fungal composition. 



Leaf Roots 



Assemblage 




Df 


F (pvalue) 


R2 


F (pvalue) 


R2 


Fungi 


region 


2 


9.768 (0.001 ) 


0.263 


4.903 (0.001 ) 


0.1 74 




pH 


1 


3.867 (0.003) 


0.052 


5.1 09 (0.001 ) 


0.091 




temperature 


1 


10.749 (0.001) 


0.145 


4.858 (0.001) 


0.086 




region:site 


7 


2.744 (0.001) 


0.258 


2.238 (0.001) 


0.278 




Residuals 


21 




0.282 




0.372 


Ascomycetes 


region 


2 


10.470 (0.001) 


0.274 


5.326 (0.001) 


0.188 




pH 


1 


4.064 (0.001) 


0.053 


5.394 (0.001) 


0.095 




temperature 


1 


1 1.001 (0.001) 


0.144 


5.787 (0.001) 


0.102 




region:site 


7 


2.765 (0.001) 


0.253 


1.994 (0.001) 


0.246 




Residuals 


21 




0.275 




0.370 


Basidiomycetes 


region 


2 


7.836 (0.001 ) 


0.228 


3.980 (0.001 ) 


0.146 




dH 


1 


3.824 (0.002) 


0.056 


4.703 (0.001) 


0.086 




temperature 


1 


10.496 (0.001) 


0.153 


3.933 (0.001) 


0.072 




region:site 


7 


2.526 (0.001) 


0.257 


2.405 (0.001) 


0.309 




Residuals 


21 




0.306 




0.386 


Ectomycorrhizal 


region 


2 






3.724 (0.001) 


0.140 


fungi 


PH 


1 






4.665 (0.001) 


0.088 




temperature 


1 






4.247 (0.001) 


0.080 




region:site 


7 






2.270 (0.001) 


0.298 




Residuals 


21 








0.394 



Permutational multivariate analysis of variance of the compositional dissimilarity. 
doi:1 0.1 371 /joumal.pone.01 00668.t004 



Furthermore, we observed that the fungal richness does not 
correlate to the soil pH (Table 3). This finding corroborates a 
similar observation of a non-significant relationship between the 
soil pH and fungal richness [24]. 

The composition of leaf-associated fungal assemblages 
covaries with temperature 

Our results confirm that mean annual temperature could be of 
main importance in structuring the phyllosphere fungal assem- 
blages at large geographic scales [1 1] (Table 4). Only one sampling 
in summertime has been done to avoid taking the seasonal 
variability into account [11]. Care should be taken about the 
conclusions as factors such as total rainfall or soil phosphorus 
content were strongly correlated to the temperature and thus, the 
exact causal factor could not be determined. Some recent studies 
showed that the leaf-associated fungal assemblages are spatially 
structured, from the regional scale [9] to the single tree canopy 
scale [48-49] and along elevation gradients [50]. Although these 
fungi are often generalists with a cosmopolitan distribution, these 
assemblages are structured by both abiotic factors such as the 
mean annual temperature [11] or rainfall [50], and biotic factors 
such as the host genotype [49]. The difference in fungal 
assemblages between, in one hand, the Pyrenees and, in another 
hand, the Alps and the Vosges could be the result of a higher initial 
sequencing effort in the Pyrenees although all the assemblages 
were randomly downsampled at a similar sequencing depth. 



The composition of root-associated fungal assemblages 
covaries both with soil pH and temperature 

Several recent studies showed that the belowground fungal 
composition varies along elevation gradients [12], [19], [20], [22]. 
In this study, the composition of root-associated fungi first correlat 
with the soil pH and secondly with the temperature (Table 4). 
However, the soil pH correlated with the C:N ratio and because 
the mean annual temperature correlated with two other soil 
variables (nitrogen and phosphorus) the direct effect of climate 
cannot be ascertained (Figure S2 and Table S3). Nevertheless, the 
soil characteristics were previously reported as drivers of the 
microbial assemblage. Indeed, our results show that fungal 
composition is strongly related to soil pH confirming previous 
reports on soil fungal and bacterial communities [24], [51], 

Our results suggest that temperature might be an important 
factor in shaping EcM specific composition [12], [16] although it 
was not possible to distinguish the effect of the climatic variables 
and the correlated soil variables. The climatic variables could drive 
EcM composition indirectly by the effect of climate on vegetation 
through root status and turnover for example. Indeed, it is known 
that a major structuring factor of EcM assemblages is the host 
family [53]. The beech-dominated stands were explicitly chosen to 
limit this biotic effect. This may explain why the region effect 
(Alps, Pyrenees, Vosges) is of less importance for explaining root- 
associated basidiomycetes assemblage diversity as EcM fungi 
closely associated with their beech host represent a large part of 
this assemblage. Some EcM fungi show a high presence 
throughout the study sites such as Cenococcum spp. Its ubiquity 
within the French beech forest has already been reported and 
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Figure 1. Relationships between fungal richness and mean annual temperature. Fungal richness from the phyllosphere (A, C, E) or 
associated with the fine-roots (B, D, F), either total richness (A, B), Ascomycetes richness (C, D) or Basidiomycete richness (E, F). The richness was 
estimated from rarefaction at a sequencing depth of 1 400 and 500 sequences per leaf and root samples, respectively. Samples corresponded to the 
Alps (closed dots), the Pyrenees (filled triangles) or the Vosges (plus sign) with different colours meaning different sites. 
doi:1 0.1 371 /journal.pone.01 00668.g001 



could result from the pooling of different cryptic species as it is now 
recognized that C. geophilum is a species complex [52], each species 
possibly having a particular niche. For instance, the MOTU 
Cenococcum spl was restricted to plots where the pH was above 5. 



Conclusions 

According to our results, the above-ground and below-ground 
fungal assemblages do not follow similar environmental drivers. 
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Figure 2. Correspondence analysis representing the distribution of the samples according to the presence/absence of fungal 
Molecular Operational Taxonomic Units (MOTUs). There were 512 phyllosphere-associated MOTUs (A), 472 root-associated MOTUs (B), and 
120 root-associated EcM MOTUs (C). Samples were from the Alps (closed dots), the Pyrenees (filled triangles) or the Vosges (plus sign) with different 
colours meaning different sites. Percentage of variance into brackets. 
doi:1 0.1 371 /journal.pone.01 00668.g002 



The phyllosphere assemblage was found to be dominated by 
ascomycetes, as has already been described [9] , [54] whereas both 
ascomycetes and basidiomycetes contributed to root-associated 
fungal assemblage in a similar proportion, even if EcM fungi are 
dominated by basidiomycetes. While the phyllosphere assemblage 
appeared to be largely related to climatic variables, the root- 
associated assemblage was related to both edaphic and climatic 
variables. To go further, it appears important to analyse the data 
at lower taxonomic levels (e.g., [24]) or taking into account the 
ecological trait differences. It is possible that the fungal taxonomic 
groups are too heterogeneous to be pooled into one assemblage. 
Considering the whole assemblage might therefore blur the 
relationship of sub-assemblage with environmental gradients and 
impede our understanding of fungal community ecology. 

Supporting Information 

Figure SI Venn diagrams with the number of non 
singleton MOTUs. MOTUs were found in the different regions 
(A), from the phyllosphere at the different elevation sites (B) and from 
the root-associated assemblages at the different elevation sites (C). 
(PPT) 

Figure S2 Principal component analysis of the environ- 
mental variables. Note that the annual precipitation and 
elevation are overlapped. 
(PDF) 

Figure S3 Relationships between fungal richness and 
elevation. Fungal richness from the phyllosphere (A, C, E) or 
associated with the fine-roots (B, D, F), either total richness (A, B), 
Ascomycetes richness (C, D) or Basidiomycete richness (E, F). The 
richness was estimated from rarefaction at a sequencing depth of 
1 400 and 500 sequences per leaf and root samples, respectively. 
Samples corresponded to the Alps (closed dots), the Pyrenees (filled 
triangles) or the Vosges (plus sign) with different colours meaning 
different sites. 
(TIF) 



Table SI List of the MOTUs with the taxonomic 
assignation and classification as ectomycorrhizal fungi 

(if applicable). 
(XLS) 

Table S2 Percentage of non-singleton MOTUs shared or 
not shared between the elevation sites of each region. 

(PPTX) 

Table S3 Pearson coefficients of correlation between 
the environmental variables. Significant correlations are 
indicated with asterisks for the following thresholds: ***<0.001< 
**<0.01<*<0.05. 
(XLS) 

Table S4 Estimated richness per plot. The richness was 
estimated from rarefaction at a sequencing depth of 1 400 and 500 
sequences per leaf and root samples, respectively. See Materials 
and Methods for the detailed estimation methods. 
(XLS) 
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